
*This code must be run within the main do file. It will not run correctly on its own. 
	
	append using "$data/US_BCCA60_2040_2050.dta"

	sort fips year
	
	replace PovCentile=PovCentile[_n-1] if fips==fips[_n-1] & PovCentile==.
	
	bysort fips: egen Hist_avg_90=mean(days90plus_dry) if year<2012
	egen all_avg=mean(days90plus_dry) if year<2012
		label var all_avg "Historical average, 1986-2011"
	
	bysort year: egen future_avg=mean(days90plus_dry) if year>2012
		label var future_avg "Annual average, 2040-2050"
		
	bysort fips: egen future2=mean(days90plus_dry) if year>2012
	
	bysort fips: egen past=mean(days90plus_dry) if year<2012
	
	g future_difference=days90plus_dry-Hist_avg_90
		label var future_difference "Difference between average # of days 90+ from 86-11 and 40-50"
	
	sort fips year
	replace mean_Div4=mean_Div4[_n-1] if fips==fips[_n-1] & mean_Div4==.
	foreach m in dayimpact90 interactimpact rawdayimpact90{
		replace `m'=`m'[_n-1] if `m'==.	
	}
	
	g net_lost_payroll=(dayimpact90+(mean_Div4*interactimpact))*future_difference
		label var net_lost_payroll "Payroll lost relative to average days above 90"
		
	g lost_payroll=(dayimpact90+(mean_Div4*interactimpact))*days90plus_dry
		label var lost_payroll "Absolute lost payroll"
		
	g raw_net_lost_payroll=(rawdayimpact90)*future_difference
		label var net_lost_payroll "Payroll lost relative to average days above 90 (raw)"
		
	g raw_lost_payroll=(rawdayimpact90)*days90plus_dry
		label var lost_payroll "Absolute lost payroll (raw)"	
	
	

	
		
		
	
	sort fips year
	
	foreach k in  PovCentile {
		replace `k'=`k'[_n-1] if fips==fips[_n-1] & `k'==.
	}
	
	foreach k in PovCentile  {
		preserve
		keep if year>2039
	
			collapse (mean) days90plus_dry net_lost_payroll raw_lost_payroll raw_net_lost_payroll lost_payroll, by(`k')
			
			
			
			
				twoway  (qfitci lost_payroll days90plus_dry) (qfit raw_lost_payroll days90plus_dry, lp(dash)) ///
						(scatter lost_payroll days90plus_dry), scheme(s1color) xscale(range(30 105)) xtitle("Avg. Annual Days above 90F, 2040-2050") ytitle("Percent Loss (Payroll)")	///
						 legend(off)
						foreach p in "$output/images" {
						
						graph export "$output/images/Fig2a.png", as(png) replace 
					}
			
		restore
	}
	
	
	foreach k in  PovCentile {
		preserve
		keep if year>2039
			collapse (mean) net_lost_payroll raw_lost_payroll raw_net_lost_payroll lost_payroll, by(`k')
			foreach i in net_lost_payroll raw_lost_payroll raw_net_lost_payroll lost_payroll{
				label var `i' "Percentage Lost Payroll"
			}
			twoway (qfitci lost_payroll `k') (qfit raw_lost_payroll `k', lp(dash))(scatter lost_payroll `k')   ///
						 , scheme(s1color) ytitle("Percent Loss (Payroll)")	///
						 legend(off)
					foreach p in $output\images{
						
						graph export "$output/images/ 	Fig2b.png", as(png) replace 
					}
			
		restore
	}
	

	//Sum stats for lost payroll
	
		preserve
		keep if year>2039
		
		bysort fips: egen avg_raw_loss=mean(raw_lost_payroll)
			label var avg_raw_loss "Annual avg. lost payroll 2040-2050 (raw)"
		
		bysort fips: egen avg_loss=mean(lost_payroll)
			label var avg_loss "Annual avg. lost payroll 2040-2050"
			
		bysort fips: egen avg_raw_net_loss=mean(raw_net_lost_payroll)
			label var avg_raw_loss "Annual avg. net lost payroll 2040-2050 (raw)"
		
		bysort fips: egen avg_net_loss=mean(net_lost_payroll)
			label var avg_net_loss "Annual avg. net lost payroll 2040-2050"	
		
		bysort fips: egen avg_days=mean(days90plus_dry)
			label var avg_days "Average days annually above 90F 2040-2050"
		
		duplicates drop fips, force
		
		g avg_raw_loss_poor=100*avg_raw_loss if PovCentile>=90
		g avg_raw_loss_rich=100*avg_raw_loss if PovCentile<=10
		g avg_loss_poor=100*avg_loss if PovCentile>=90
		g avg_loss_rich=100*avg_loss if PovCentile<=10
			
		g avg_raw_net_loss_poor=100*avg_raw_net_loss if PovCentile>=90
		g avg_raw_net_loss_rich=100*avg_raw_net_loss if PovCentile<=10
		g avg_net_loss_poor=100*avg_net_loss if PovCentile>=90
		g avg_net_loss_rich=100*avg_net_loss if PovCentile<=10
		
		g avg_days_poor=avg_days if PovCentile>=90
		g avg_days_rich=avg_days if PovCentile<=10
		
		replace avg_loss_rich=11.18249 if avg_loss_rich>11.183 & avg_loss_rich!=.
		replace avg_loss_poor=-11.17784 if avg_loss_poor<-11.17784 & avg_loss_poor!=.
		
		
				rename (avg_days_poor avg_days_rich) (Poor Rich)
					label var Poor "Poor Counties"
					label var Rich "Wealthy Counties"
					cap estpost su Poor Rich
					est store days
					mat c=e(mean)
					local zpoor=c[1,1]
					local zrich=c[1,2]
				rename (Poor Rich) (avg_days_poor avg_days_rich)
			
			rename (avg_raw_loss_poor avg_raw_loss_rich) (Poor Rich)
				label var Poor "Poor Counties"
				label var Rich "Wealthy Counties"
				cap estpost su Poor Rich
				est store naiveloss
				mat a=e(mean)
					local xpoor=a[1,1]
					local xrich=a[1,2]
			rename (Poor Rich) (avg_raw_loss_poor avg_raw_loss_rich)
			
			
			rename (avg_loss_poor avg_loss_rich) (Poor Rich)
				label var Poor "Poor Counties"
				label var Rich "Wealthy Counties"
				cap estpost su Poor Rich
				est store loss
				mat b=e(mean)
					local ypoor=b[1,1]
					local yrich=b[1,2]
			rename (Poor Rich) (avg_loss_poor avg_loss_rich)
			
			local ratiopoor=`ypoor'/`xpoor'
			local ratiorich=`yrich'/`xrich'
			
			g ratiopo=`ratiopoor'
			g ratiori=`ratiorich'
			
			rename (ratiopo ratiori) (Poor Rich)
				label var Poor "Poor Counties"
				label var Rich "Wealthy Counties"
				cap estpost su Poor Rich
				est store ratio
			rename (Poor Rich) (ratiopo ratiori)
			
			g daysdiff=`zpoor'-`zrich'
			g naivediff=`xpoor'-`xrich'
			g adjdiff=`ypoor'-`yrich'
			
			rename (daysdiff) (diff)
				label var diff "Difference"
				estpost su diff
				est store diff1
			rename (diff) (daysdiff)
			
			rename (naivediff) (diff)
				label var diff "Difference"
				estpost su diff
				est store diff2
			rename (diff) (naivediff)
			
			rename (adjdiff) (diff)
				label var diff "Difference"
				estpost su diff
				est store diff3
			rename (diff) (adjdiff)
			
			
		

	
		
		
		esttab days naiveloss loss ratio using "$output/table_projection.tex", replace ///
			refcat(label) ///
			mtitle("Days over 90$\degree\$F"  "Na{i}ve (\%)" "Adjusted Loss (\%)" "\shortstack{Ratio\\(Adj. Loss/Na{i}ve Loss)}") ///
			cells(mean(fmt(1))) label booktabs nonum collabels(none) f noobs
			
		esttab diff1 diff2 diff3 using "$output/table_projection.tex", append ///
			refcat(label) ///
			cells(mean(fmt(1))) label booktabs nonum collabels(none) f noobs
		
		restore
